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ABSTRACT 

We  improve  the  performance  of  the  Perfectly  Matched  Layer  by  using  an  automatic  fep-adaptive  discretization. 
By  means  of  /ip-adaptivity,  we  obtain  a  sequence  of  discrete  solutions  that  converges  exponentially  to  the 
continuum  solution.  Asymptotically,  we  thus  recover  the  property  of  the  PML  of  having  a  zero  reflection 
coefficient  for  all  angles  of  incidence  and  all  frequencies  on  the  continuum  level.  This  allows  us  to  minimize 
the  reflections  from  the  discrete  PML  to  an  arbitrary  level  of  accuracy  while  retaining  optimal  computational 
efficiency.  Since  our  Zip-adaptive  scheme  is  automatic,  no  interaction  with  the  user  is  required.  This  renders 
tedious  parameter  tuning  of  the  PML  obsolete.  We  demonstrate  the  improvement  of  the  PML  performance 
by  Zip-adaptivity  through  numerical  results  for  acoustic,  elastodynamic  and  electromagnetic  wave-propagation 
problems  in  the  frequency  domain  and  in  different  systems  of  coordinates. 

Keywords  and  Phrases:  Perfectly  Matched  Layer,  Zip-adaptivity,  exterior  boundary-value  problem,  acoustic 
scattering,  elastodynamic  wave  propagation,  electromagnetic  scattering. 


1.  Introduction 

Most  wave-propagation  problems  arising  in  acoustics,  elastodynamics  and  electromagnetics  are  posed 
on  a  spatially  unbounded  domain.  Since  domain-based  discretization  methods  such  as  finite  elements 
can  handle  only  bounded  domains,  a  truncation  of  the  unbounded  physical  domain  to  a  bounded 
computational  domain  is  necessary.  Any  artificial  boundary  condition  that  is  applied  on  the  truncated 
domain  boundary  must  not  essentially  alter  the  original  problem,  i.e.  it  must  allow  for  outgoing 
waves  only  and  it  has  to  be  essentially  reflectionless.  Many  techniques  have  been  developed  for  this 
purpose  such  as  Infinite  Elements  (see  e.g.  [4,  12]),  the  Dirichlet-to-Neumann  map  (see  e.g.  [14]), 
exact  nonreflecting  boundary  conditions  (see  e.g.  [15]),  the  continued-fraction  absorbing  boundary 
condition  [16]  and  the  Perfectly  Matched  Layer  [3].  We  refer  to  Refs.  [25,  30]  for  an  overview  over 
these  so-called  absorbing  boundary  conditions,  all  of  which  have  their  specific  strengths  and  weaknesses. 

Among  the  various  absorbing  boundary  conditions,  the  Perfectly  Matched  Layer  (PML)  in  particu¬ 
lar  has  become  very  popular  due  to  its  performance,  conceptual  simplicity  and  ease  of  implementation. 
The  PML  is  based  on  the  concept  of  analytic  continuation  of  a  real  function  into  the  complex  plane 
(see  e.g.  [13]),  typically  also  referred  to  as  complex- coordinate  stretching  [8].  This  concept  enables 
straightforward  derivation  of  the  PML  equations  for  different  types  of  wave  propagation  and  differ¬ 
ent  systems  of  coordinates;  see  e.g.  Refs.  [23,  24,  29,  31].  The  PML  concept  has  been  extensively 
studied  and  analyzed  regarding  the  so-called  split  and  unsplit  formulations  and  well-posedness  (see 
e.g.  Refs.  [1,  18]),  causality  [6,  28],  long-time  behavior  [1,  2],  termination  of  the  PML  [9,  26]  and 
finite-element  dispersion  and  convergence  analysis  in  the  frequency  domain  [5,  17].  A  recommendable 
review  of  the  PML  in  the  context  of  Maxwell’s  equations  is  presented  in  Ref.  [29]. 

The  PML  has  the  remarkable  property  of  having  a  zero  reflection  coefficient  for  all  angles  of 
incidence  and  all  frequencies  on  the  continuum  level.  However,  under  discretization,  this  property  is 
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compromised  and  spurious  reflections  typically  do  occur  at  the  discrete  PML.  Much  effort  has  thus 
been  spent  on  optimizing  the  PML  parameters  to  minimize  the  error  resulting  from  discretization 
and  thus  obtain  optimal  performance  of  the  PML;  see  e.g.  [7,  9,  27].  In  this  work,  we  advocate 
the  opposite  approach,  i.e.  we  optimize  the  discretization  for  any  given  PML  parameters,  rather 
than  resorting  to  tedious  parameter  tuning.  A  highly  effective  methodology  to  construct  an  optimal 
discretization  is  the  /ip-adaptive  finite-element  method  that  allows  for  locally  adapting  the  element 
size  h  and  polynomial  approximation  order  p  to  the  local  resolution  requirements  of  the  solution.  Such 
adaptivity  can  be  accomplished  in  an  automatic  way,  i.e.  no  interaction  with  the  user  is  required. 
This  is  achieved  by  a  so-called  two-grid  strategy  [10]  that  guides  the  adaptive  refinements.  By  solving 
the  problem  on  a  sequence  of  coarse  and  corresponding  fine  grids,  an  approximation  of  the  error 
function  is  constructed  at  each  iteration,  which  serves  as  an  error  estimate  and  based  on  which  it 
is  determined  where  to  refine  and  how  to  refine,  i.e.  in  h  or  in  p.  This  yields  optimal  accuracy  for 
the  least  computational  expense.  For  details  on  the  automatic  ft-p-adaptive  procedure,  we  refer  to 
Ref.  [10]  and  to  the  recent  book  by  Demkowicz  [11].  The  immense  benefit  that  ft-p-adaptivity  can 
offer  for  increasing  the  performance  of  the  PML  is  due  to  the  property  of  the  PML  of  having  a  zero 
reflection  coefficient  on  the  continuum  level.  Although  this  property  is  lost  under  discretization,  it  can 
be  recovered  asymptotically,  i.e.  upon  convergence  to  the  continuum  solution.  This  implies  that  upon 
convergence  of  the  /ip-adaptive  discretization,  the  discretization  error  and  the  resulting  reflections  can 
be  minimized  to  an  arbitrary  level  of  accuracy. 

The  PML  is  designed  such  that  traveling  waves  that  leave  the  actual  domain  of  interest  and  enter 
the  encompassing  PML  are  converted  into  evanescent  waves  that  are  made  to  decay  sufficiently  fast, 

1. e.  they  vanish  before  they  can  get  reflected  at  the  domain  boundary.  The  rapid  decay  that  is  enforced 
on  the  solution  in  the  PML  implies  strong  solution  gradients  which  need  to  be  resolved  to  prevent 
reflections  and  an  accompanying  loss  in  accuracy.  Although  the  significance  of  resolving  the  PML  has 
long  been  recognized,  conventional  discretization  methods  generally  face  difficulties  in  resolving  strong 
solution  gradients  and,  therefore,  typically  revert  to  adjusting  the  PML  parameters  and  the  associated 
damping  profile  to  a  given  discretization.  Conversely,  our  adaptive  strategy  automatically  adapts  the 
discretization  to  resolve  any  damping  profile  that  is  imposed  on  the  solution.  This  allows  us  to  choose 
the  geometry  of  the  PML  independent  of  the  geometry  of  the  initial  grid,  i.e.  the  respective  geometries 
need  not  coincide  (a  cylindrical  PML  and  a  Cartesian  grid,  for  instance)  -  a  truly  remarkable  feature.  In 
this  paper,  we  show  that  by  means  of  /ip-adaptivity  PML-induced  solution  gradients  can  be  effectively 
resolved  and  the  performance  of  the  PML  can  be  greatly  improved.  This  renders  tedious  parameter 
tuning  of  the  PML  obsolete.  We  demonstrate  the  performance  and  the  versatility  of  our  approach  by 
numerical  results  for  acoustic,  elastodynamic  and  electromagnetic  wave-propagation  problems  in  the 
frequency  domain  and  in  different  systems  of  coordinates. 

The  contents  of  this  paper  are  organized  as  follows:  Sections  2,3  and  4  present  the  PML  formulation 
of  the  Helmholtz,  linear  elasticity  and  Maxwell’s  equations,  respectively,  and  discuss  the  particularities 
of  the  respective  PML  formulation.  Each  of  these  sections  presents  numerical  results  that  demonstrate 
the  potential  of  /ip-adaptivity  for  improving  the  performance  of  the  PML.  Finally,  Section  5  contains 
concluding  remarks. 

2.  PML  FORMULATION  FOR  THE  HeLMHOLTZ  EQUATION 

In  this  section,  we  study  the  formulation  of  the  Perfectly  Matched  Layer  for  an  exterior  boundary- 
value  problem  governed  by  the  Helmholtz  equation.  In  Section  2.1,  we  provide  a  problem  statement 
in  coordinate-invariant  form.  In  Section  2.2,  we  derive  the  PML  formulation  in  Cartesian  coordinates. 
Based  on  the  formalism  introduced  in  Section  2.2,  in  Sections  2.3  and  2.4,  we  state  the  PML  formu¬ 
lation  in  2D  polar  and  3D  spherical  coordinates,  respectively.  In  Section  2.5,  we  present  numerical 
results  for  exterior  Helmholtz  problems  that  demonstrate  the  improvement  in  performance  of  the  PML 
that  can  be  obtained  by  /ip-adaptive  finite  elements. 
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2.1  Problem  statement 

We  consider  the  unbounded  exterior  domain  Q,  =  M"— fi'"*  with  fi'”*  C  M",  n  =  1,  2, 3,  a  given  bounded 
interior  domain  and  n  the  number  of  space  dimensions.  We  seek  a  solution  for  the  pressure  p{x),  a;  G 
that  satisfies  the  Helmholtz  equation  in 

—Ap  —  k^p  =  0  infl,  (2.1) 


where  A  denotes  the  Laplacian  operator,  and  k  :=  to  jc  is  the  wavenumber  with  oj  and  c  the  angular 
frequency  and  the  sound  speed,  respectively.  Throughout  this  work,  we  shall  assume  a  time-harmonic 
dependence  of  the  form  where  i  is  the  imaginary  unit  and  t  denotes  time  and,  moreover,  we 
consider  the  governing  equations  in  their  non-dimensional  form.  Eq.  (2.1)  is  supplemented  with  a 
Dirichlet  boundary  condition  on  the  interior-domain  boundary 

p  =  Pd  onr  =  9n‘”*^,  (2.2) 

and  the  Sommerfeld  radiation  condition  at  infinity 

^  +  ikpeL^{n),  (2.3) 


where  r  denotes  the  radius.  Note  that  this  particular  form  of  the  Sommerfeld  radiation  condition  has 
the  advantage  of  being  independent  of  the  number  of  space  dimensions.  We  further  remark  that  the 
Dirichlet  boundary  condition  (2.2)  can  be  replaced  by  a  Neumann  or  a  Cauchy  (impedance)  boundary 
condition. 


dp 

fW  = 

on 


(2.4a) 


dp 

dn 


+  iujf3p  =  g  , 


(2.4b) 


on  the  entire  boundary  or  on  parts  of  it,  where  n  denotes  the  outward  unit  normal  on  T,  and  (3  >  Q 
and  g  are  given  functions  on  T. 

Finally,  let  us  remark  that  the  Helmholtz  equation  is  obtained  by  eliminating  the  velocity  from 
the  first-order  system  of  linear-acoustics  equations  which,  in  coordinate-invariant  form,  are 


iutp  +  poc^ V  ■  u  =0 
iiopou  +  Vp  =  0 , 


(2.5) 


where  p  and  u  denote  the  small-amplitude  disturbance  in  pressure  and  velocity,  respectively,  of  a 
uniform  stationary  fluid,  and  po  is  the  fluid  density  associated  with  this  uniform  state. 


2.2  PML  formulation  in  Cartesian  coordinates 

The  governing  equations  in  the  PML  can  be  derived  by  applying  a  complex  coordinate  transformation 
of  the  original  equations  in  the  direction  normal  to  the  interface  between  the  “domain  of  interest” 
and  the  PML.  For  this  purpose,  it  is  essential  that  the  solution  to  the  considered  equations  is  analytic 
in  the  direction  normal  to  the  interface.  Analyticity  of  the  solution  is  equivalent  to  the  solution 
being  holomorphic,  i.e.  differentiable  in  the  complex  sense.  This  assumption  allows  us  to  consider 
the  solution  as  a  function  of  complex  variables  Zi  in  place  of  real  variables  Xi,  and  to  replace  the 
derivative  with  respect  to  Xi  with  the  derivative  with  respect  to  Zi.  This  procedure  is  commonly  also 
referred  to  as  analytic  continuation  [13]  or  as  complex- coordinate  stretching  [8].  The  motivation  of 
such  coordinate  stretching  is  to  construct  a  continuation  of  the  solution  into  the  complex  plane  in 
such  a  way  that,  in  the  PML,  the  solution  becomes  evanescent,  whereas  in  the  domain  of  interest  the 
original  solution  is  retained. 
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The  analytic  continuation  of  the  solution  into  the  complex  plane  is  obtained  by  the  following 
transformation  of  the  j-th  coordinate 


■■=  [xj  +  aj{xj)]  -  i[bj{xj)]  j  = 


(2.6a) 


d  d  Id  .  ,  /  dzjixj) 

T. - >  t  —  ^  ^3  ■= 

axj  azj{xj)  Zj  oxj  ■’  dxj 


1— ^ 

■f 

—  i 

'dbjixjY 

L  dxj  \ 

dxj 

(2.6b) 


where  the  j-th  complex  coordinate  Zj  depends  only  on  the  j-th  real  coordinate  Xj.  Factors  aj{xj)  and 
bj{xj)  are  “suitable”  functions  specified  in  the  sequel.  In  the  domain  of  interest,  aj{xj)  =  bj{xj)  =  0 
which  yields  Zj{xj)  =  Xj  and,  thus,  the  original,  “unstretched”  equations  are  recovered.  A  jump  in  Oj, 
bj  across  the  interface  between  the  domain  of  interest  and  the  PML  induces  a  loss  of  regularity.  To 
minimize  such  loss  of  regularity,  aj(xj)  and  bj{xj)  are  typically  set  to  increase  steadily  from  zero  as 
one  moves  away  from  the  interface  through  the  PML,  ideally  with  C^-continuity  across  the  interface. 
Considering  a  wave  traveling  to  the  right,  in  the  PML,  aj(xj)  >  0  to  ensure  that  evanescent  waves 
have  an  exponential  decay  that  is  faster  in  the  PML  than  in  the  domain  of  interest;  and  bj{xj)  >  0  to 
ensure  that  traveling  waves,  once  they  enter  the  PML,  are  converted  into  evanescent  waves  that  decay 
exponentially.  Conversely,  for  a  wave  traveling  to  the  left,  we  need  aj{xj)  <  0  and  bj(xj)  <  0  in  the 
PML.  We  remark  that,  in  general,  bj{xj)  is  chosen  frequency-dependent;  see  for  instance  Ref.  [29]. 
However,  in  this  work,  we  shall  not  consider  frequency-dependent  stretching,  since  all  our  examples 
consider  a  single  frequency  only.  By  an  appropriate  choice  of  aj{xj)  and  bj{xj),  a  sufficiently  fast 
decay  of  the  solution  can  be  enforced  such  that  the  solution  of  the  PML-modified  problem  practically 
vanishes  at  the  truncated  domain  boundary  and,  thus,  satisfies  a  homogeneous  Dirichlet  boundary 
condition.  Note,  however,  that  the  enforced  rapid  decay  of  the  solution  in  the  PML  induces  strong 
solution  gradients.  We  demonstrate  in  Section  2.5  the  importance  of  resolving  these  PML-induced 
solution  gradients  for  the  accuracy  of  the  solution  in  the  domain  of  interest  and,  moreover,  that  this 
can  be  efficiently  achieved  by  hp-adaptivity.  Finally,  let  us  remark  that  the  choice  of  aj{xj)  and  bj{xj) 
has  to  respect  the  principle  of  causality;  see  Refs.  [6,  28]  for  details. 

In  the  sequel,  we  derive  the  PML  formulation  for  the  Helmholtz  problem  considered  in  Section  2.1 
in  n  dimensions  and  Cartesian  coordinates.  To  this  end,  we  start  from  the  time-harmonic  linear- 
acoustics  equations  (2.5)  expressed  in  Cartesian  coordinates 


n  r. 

9  UUj 

+  Poc  =  0 , 

i=i  ^ 

,  dp 

tupoUj  +  - —  =  0 . 


(2.7) 


Introducing  in  Eqs.  (2.7)  complex-coordinate  stretching  according  to  (2.6)  with  Xj  Zj{xj),  we  obtain 


ILOp  - 


2  1  duj 


j=l  ^3  ^^3 
1 


lUipoUj 


dp 
z';  dXi 


=  0, 

=  0. 


(2.8) 


Note  that  we  do  not  use  the  Einstein  summation  convention  here.  To  express  Eq.  (2.8) i  in  a  weak 
form,  we  multiply  (2.8) i  with  z'  :=  and  then  with  a  test  function  q  €  Q  and  integrate  over 

the  domain  f2.  Thus,  we  obtain 


r  r  z'  du- 

iijj  /  z'pqdVl  +  pq(?22,  /  ~Ta~^dd^  =  0.  (2.9) 

Jfl  j=ld^^3d^3 

Note  that  in  deriving  the  weak  form  (2.9)  the  multiplication  of  the  expression  by  z'  facilitates  the 
subsequent  integration-by-parts,  because  the  prefactor  z' jz'^  in  (2.9)  is  independent  of  Xj.  Upon 
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integrating  Eq.  (2.9)  by  parts,  we  multiply  the  resulting  expression  by  iuj  and  eliminate  Uj  by  invoking 
Eq.  (2.8)2.  This  yields  the  final  variational  form 


f  pGpD  +  Q 


< 


z'  dp  dq 
{z'r  dx,  dx, 


Vge  Q, 


(2.10) 


where  k  =  uj/c  is  the  wavenumber  and  g  denotes  the  normal  pressure  gradient  specified  on  the 
Neumann  portion  of  the  boundary,  E jv-  In  (2.10),  pu  denotes  a  finite-energy  lift  of  the  Dirichlet  data, 
Q  is  the  space  of  test  functions 


Q-.=  {qeX  :  g  =  0onr£,}. 


and  the  “energy  space”  X  associated  with  the  variational  problem  (2.10)  is  defined  as 


T  := 


dq 

dxj 


G 


(2.11) 


(2.12) 


Note  that  in  deriving  Eq.  (2.10),  Eq.  (2.7)i  has  been  treated  in  its  weak  form,  whereas  Eq.  (2.7)2  has 
been  treated  in  its  strong  form.  This  has  implications  for  the  regularity  requirements  of  the  solution 
variables.  In  particular,  this  relaxes  the  regularity  requirements  for  the  velocity  Uj. 

Finally,  let  us  remark  that,  since  the  complex-coordinate  stretching  according  to  (2.6)  is  straightfor¬ 
ward  for  first-order  derivatives,  the  variational  formulation  of  the  PML  equations  is  most  conveniently 
derived  by  starting  from  a  system  of  first-order  equations  and  following  the  procedure  described  above. 
This  equally  applies  to  the  PML  formulation  of  the  linear-elasticity  equations  as  well  as  to  the  PML 
formulation  of  Maxwell’s  equations  that  are  stated  in  Sections  3.1  and  4.1,  respectively. 


2.3  2D  Helmholtz  equation  in  polar  coordinates 

In  this  section,  we  present  the  PML  formulation  for  the  Helmholtz  equation  in  two  dimensions  and 
polar  coordinates  {r,9)  with  x  =  rcosO  and  y  =  rsinO.  We  employ  complex-coordinate  stretching 
only  in  the  r  coordinate,  i.e.  normal  to  the  interface  between  the  domain  of  interest  and  the  PML. 
Thus,  the  solution  to  the  Helmholtz  equation  is  required  to  be  analytic  in  the  r  coordinate  only,  but 
not  necessarily  in  the  0  coordinate. 

Recalling  the  formulas  for  the  gradient  of  a  scalar  p  and  for  the  divergence  of  a  vector  field 
u  =  UrBr  +  ugEg  in  polar  coordinates 


_  dp  1  dp 


(2.13a) 


19  1  dug 

V  •  M  =  -  —  {rur)  H - ^  , 

r  dr  r  dd 


(2.13b) 

respectively,  we  convert  the  first-order  system  of  linear-acoustics  equations  (2.5)  to  polar  coordinates 


iojp  +  poc 


Id  I dug\ 

,  dp  n 

ILUpoUr  +7^  =  0 

or 

,  1 

lUjpoUg  +  -  7^  =  0  . 

r  oO 


(2.14) 
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Introducing  in  (2.14)  complex-coordinate  stretching  according  to  (2.6)  with  r  z{r)  yields 


iujp  +  poc 


Id.  ,  1  du0 

—  —  {zUr)  H - ^ 

z'z  or  z  da 


\z'zdr''~^"  z  d9  J 

<  iujpoUr+^^  =0  (2-15) 

z  Or 

dp 

lUJpQUg  +  -  ^  =  0  . 

z  do 

Multiplying  Eq.  (2.15)i  by  iwz' zjr  and  subsequently  by  a  test  function  q  G  Q,  integrating  the  resulting 
expression  over  domain  and  integrating  by  parts,  we  obtain 


iLopoug 


-pq  rdrdO  —  iujpoc^ 


f  { z  dq  z'  dq 


rdrdO  =  0 , 


where  we  tacitly  assumed  Dirichlet  boundary  conditions  for  the  solution  on  the  entire  boundary,  E, 
and,  accordingly,  homogeneous  boundary  conditions  for  the  test  function  q.  To  eliminate  the  velocity 
components  Ur  and  ug  from  Eq.  (2.16),  we  substitute  Eqs.  (2.15)2,3  into  (2.16),  and  obtain  the  final 
variational  formulation 


p  G  Pd  +  Q 

f  + -II) 

Jn  \  ^  f  dr  Or  rz  00  00  J 


/  z  z 

/  - pq  rdrdO  =  0 

7  ^ 


I  Vg  G  Q . 

In  (2.17),  pd  denotes  a  finite-energy  lift  of  the  Dirichlet  data,  Q  is  the  space  of  test  functions 

Q:={qGX  :  g  =  0onr}, 

and  the  “energy  space”  X  associated  with  the  variational  problem  (2.17)  is  defined  as 


Since  the  solution  to  problem  (2.17)  decays  exponentially  in  r,  we  replace,  upon  discretization,  the 
infinite  physical  domain  by  a  finite  computational  domain.  We  can  then  discretize  problem  (2.17)  on 
the  truncated  domain  by  means  of  finite  elements.  Note  that  the  truncated  domain  need  not  have  a 
circular  shape;  moreover,  we  typically  revert  to  Cartesian  coordinates  Xj  to  represent  the  integrand 
of  Eq.  (2.17).  In  particular,  the  integrand  of  the  first  integral  in  (2.17)  can  then  be  rewritten  as 

z  dp  dq  z'  dp  dq  z  dpdq  z'r  /  dp  dq  dp  dq\ 

z'r  dr  dr  rz  dO  dO  z'r  dr  dr  z  \  dxj  dxj  dr  dr ) 

z'r  dp  dq  (  z  z'r\  dpdq  /o 


z  dxj  dxj 
dp  dq 
dxj  dxi 


(2.20a) 


z  \  z'r  z 


z  r  \  XiXi 


(2.20b) 


where  6ij  denotes  the  Kronecker  delta  and  the  Einstein  summation  convention  has  been  used.  The 
first  equality  in  (2.20a)  is  obtained  by  invoking  Eq.  (2.13a),  and  the  third  equality  is  obtained  by 
using  the  fact  that 

dr  dx,  dr  dx,  r  '  ^  ^ 


and  likewise  for  dq/dr. 
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2.^  3D  Helmholtz  equation  in  spherical  coordinates 

In  this  section,  we  present  the  PML  formulation  for  the  Helmholtz  equation  in  three  dimensions  and 
standard  spherical  coordinates  (r,  with  x  =  rsmipcosO,  y  =  rsin^/;sin0  and  2;  =  r  cos  ip.  The 
derivation  proceeds  analogously  to  the  case  of  polar  coordinates;  cf.  Section  2.3. 

Recalling  the  formulas  for  the  gradient  of  a  scalar  p  and  the  divergence  of  a  vector  field  u  = 
UrEr  +  +  ugeg  in  spherical  coordinates 


—  dp  I  dp  I  dp 

dr  r  dtp  r  sm  tp  df) 


_  ^  ^  r  2  \ 

V  ■u  =  -^^(r  u^)  + 
or 


1  d  -It  .  1  9u0 

r  sm  tp  dtp  r  sm  tp  da 


(2.22a) 

(2.22b) 


we  convert  the  first-order  system  of  linear-acoustics  equations  (2.5)  to  spherical  coordinates 


2  i  ^  d  2  \  1  ^ 

lojp  +  poc  \  Ur)^ - ^ 

or  r  sm  ip  oip 


{up,  sin  tp)  ■ 


ZUJ  Pq  Ujj' 

iujpoup,  + 
iojpoue 


1  dug 
r  sin  tp  do 

dp 

1  dp 
r  dtp 
1  dp 
r  sin  tp  do 


=  0 
=  0 
=  0 
=  0. 


(2.23) 


In  spherical  coordinates,  the  reflectionless  absorption  of  outward  traveling  waves  can  be  achieved 
through  analytic  continuation  on  the  radial  variable  r  according  to  Eq.  (2.6)  with  r  ^  z{r).  This 
transforms  Eq.  (2.23)  into 


iojp  +  poc 


1  i9  2  ^  1  '9 

-(Z  Ur)  + 


'  z^  dr 


z  sin  tp  dtp 


{up,  sin  tp) 


1  dug 
z  sin  tp  do 

1  dp 

lUtpoUr  + 

z'  dr 
1  dp 
z  dtp 
1  dp 
z  sin  tp  do 


lUjpoUp, 

iutpoug 


=  0 
=  0 
=  0 
=  0. 


(2.24) 


We  multiply  Eq.  (2.24)  1  by  iioz' z^  jr"^  and  then  by  a  test  function  q  G  Q,  and  we  integrate  the  resulting 
expression  over  the  domain  f2,  which  yields 


—  ui 


z'  -^pq  sin  tpdrdtpdO 


■  lujpoc 


m  r‘ 

[ 


In 


15/2  \ 

—  —  (z  Ur)q  + 
r^  or 


d 


sin  tp  dtp 


{up,  sin  tp)q  - 


dug 


sin  tp  do 


q  sin  tpdrdtpdO  =  0  .  (2.25) 


Upon  integrating  Eq.  (2.25)  by  parts  and  eliminating  Ur,  up,  and  ug  by  invoking  Eqs.  (2.24)2,3,4, 
respectively,  we  obtain  the  final  variational  formulation 


P&Pd  +  Q 


f  P  z"^  dp  dq  z'  dp  dq 

Jn  xz'r"^  dr  dr  dtp  dtp 


z'  dp  dq 
sin^  tp  do  dO 


sin  tpdrdtpdO 


Vge  Q, 


(2.26) 
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where  g  denotes  the  normal  pressure  gradient  specified  on  the  Neumann  portion  of  the  boundary,  Fjv. 
For  the  remainder  of  the  boundary,  T]j,  we  assume  homogeneous  Dirichlet  boundary  conditions  for 
the  pressure  and,  accordingly,  homogeneous  boundary  conditions  for  the  test  function  q.  In  (2.26), 
Pd  denotes  a  finite-energy  lift  of  the  Dirichlet  data,  Q  is  the  space  of  test  functions 

Q:={qeX  :  g  =  0  on  F^}  ,  (2.27) 


and  the  “energy  space”  X  is  defined  as 


z'\2r  or  r 


where  Vg  denotes  the  gradient  on  the  unit  sphere 


_  I  dq  dq 

sm  oO  dip 


(2.28) 


(2.29) 


2.5  Numerical  experiments 

To  demonstrate  the  benefit  of  /ip-adaptivity  for  the  performance  of  the  PML,  we  present  below  nu¬ 
merical  results  of  wave  propagation  in  an  acoustic  medium.  In  particular,  we  consider  the  scattering 
of  a  2D  plane  wave  on  a  cylinder  with  a  polar  PML  and,  subsequently,  the  scattering  of  a  3D  plane 
wave  on  a  sphere  with  spherical  PML. 


2.6  2D  polar  PML 

In  this  section,  we  test  the  /ip-adaptive  algorithm  for  a  two-dimensional  PML  in  polar  coordinates 
for  the  time-harmonic  linear-acoustics  equations  according  to  Eq.  (2.17).  This  problem  is  suitable  for 
assessing  the  performance  of  the  PML  combined  with  /ip-adaptivity  and  the  accuracy  of  the  numerical 
solution,  since  the  problem  admits  an  analytic  solution  (see  e.g.  [20,  ch.  10.5]),  which  is  expressed  in 
terms  of  standard  harmonics  and  Hankel  functions.  The  validity  of  the  analytic  solution  is  extended 
to  the  PML  by  means  of  complex-coordinate  stretching.  Note  that  the  computational  domain  need 
not  have  a  cylindrical  shape.  In  fact,  to  challenge  our  /ip-adaptive  algorithm,  we  select  a  cylindrical 
geometry  of  the  PML  that  does  not  match  the  Cartesian  structure  of  the  initial  mesh.  Thus,  we 
“truncate”  the  Cartesian  computational  domain  (—5,5)^  by  a  cylindrical  PML  for  r  >  4.  Complex- 
coordinate  stretching  is  invoked  in  the  r  coordinate  only,  in  particular 


z(r) 


r  0  <  r  <  4 


(2.30) 


Let  us  however  remark  that  the  stretching  (2.30)  is  just  one  possible  choice  that  ensures  a  sufficiently 
rapid  decay  of  the  solution  and  that  other  forms  are  conceivable. 

Below,  we  present  numerical  results  for  the  scattering  of  an  incident  plane  wave  on  a 

cylinder.  We  set  the  dimensionless  cylinder  radius  r  =  1,  the  dimensionless  wavenumber  k  =  Xk jZ  and 
the  direction  of  the  incident  wave  v  =  (1,0).  We  augment  the  problem  with  homogeneous  Dirichlet 
boundary  conditions  p  =  0  on  the  outer  domain  boundary  adjacent  to  the  PML  and,  moreover,  with 
Dirichlet  boundary  conditions  on  the  cylinder  boundary,  i.e.  at  r  =  1,  that  are  specified  according  to 
the  analytic  solution  from  Ref.  [20,  ch.  10.5]. 

Starting  from  the  initial  mesh  shown  in  Fig.  1  (left),  we  plot  in  Fig.  4  (right)  the  convergence 
curve  in  terms  of  percent  relative  error  in  the  iL^-seminorm  (in  a  logarithmic  scale)  versus  the  total 
number  of  degrees-of-freedom  N  (in  the  algebraic  scale  The  error  is  evaluated  over  the  entire 

computational  domain  including  the  PML.  Fig.  4  (right)  displays  two  convergence  curves  for  the 
cylindrical  PML  that  use  as  a  reference  solution  the  exact  analytic  solution  and  the  fine-grid  solution, 
respectively.  These  two  convergence  curves  are  almost  indistinguishable,  which  indicates  that  the 
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error  estimate  obtained  with  the  fine-grid  solution  as  reference  solution  is  indeed  very  accurate.  The 
/ip-refined  mesh  corresponding  to  an  error  level  of  1%  is  obtained  after  16  iterations;  see  Fig.  1  (right). 
Note  the  refinements  that  have  been  selected  by  the  adaptive  algorithm  in  particular  in  the  PML 
to  resolve  the  cylindrical  geometry  of  the  PML  and  the  solution  that  is  made  to  decay  rapidly  in 
radial  direction  in  the  PML.  Figs.  2  (left)  and  (right)  show  the  real  component  of  the  solution  and 
of  the  error  function,  respectively.  Here,  the  error  function  is  computed  using  the  analytic  solution 
to  the  problem;  cf.  [20].  These  plots  undermine  that  hp-adaptivity  is  capable  of  effectively  resolving 
the  rapidly  decaying  solution  in  the  PML,  thereby  minimizing  the  reflections  from  the  PML  to  an 
arbitrary  level  of  accuracy  (here  1%)  without  the  necessity  of  PML-parameter  tuning.  These  results 
are  truly  remarkable  given  the  mismatch  between  the  cylindrical  PML  and  the  Cartesian  structure  of 
the  initial  mesh. 

Clearly,  if  the  geometry  of  the  PML  is  chosen  according  to  the  structure  of  the  underlying  mesh, 
significantly  fewer  degrees-of-freedom  are  required  to  deliver  the  same  accuracy.  To  undermine  this 
claim,  we  recompute  the  above  test  case  with  a  Cartesian  PML  specified  according  to 


Zk{Xk) 


Xk  0  <  \xk\  <  4 

-  i  Xk  \xk\  >  4 


(2.31) 


Starting  from  the  same  initial  mesh  as  displayed  in  Fig.  1  (left),  an  /ip-mesh  corresponding  to  an  error 
level  of  1%  is  obtained  after  15  iterations;  see  Fig.  3  (left)  for  the  hp-refined  mesh.  Fig.  4  (left)  for  the 
real  component  of  the  error  function  with  the  fine-grid  solution  used  as  a  reference,  and  Fig.  4  (right) 
for  the  convergence  curve.  The  real  component  of  the  solution  obtained  on  the  /ip-refined  mesh  is 
shown  in  Fig.  3  (right).  Comparing  the  final  hp-meshes  for  a  cylindrical  and  a  Cartesian  PML,  Fig.  1 
(right)  and  Fig.  3  (left),  respectively,  we  observe  that  the  hp-mesh  for  a  cylindrical  PML  consists  of 
8325  degrees-of-freedom,  whereas  the  hp-mesh  for  a  Cartesian  PML  consists  only  of  4779  degrees-of- 
freedom;  compare  also  the  respective  convergence  curves  in  Fig.  4  (right).  Hence,  the  efficiency  of 
the  discretization  can  be  significantly  improved  by  constructing  initial  meshes  in  accordance  with  the 
geometry  of  the  PML. 


Figure  1:  Acoustic  scattering  of  a  plane  wave  on  a  unit  cylinder  with  a  cylindrical  PML.  Left:  Initial 
finite  element  mesh.  Right:  hp-refined  mesh  corresponding  to  1%  error  and  consisting  of  8325  degrees- 
of-freedom.  The  color  bar  indicates  the  order  of  element  edges  and  interiors. 


3D  spherical  PML 

In  this  section,  we  consider  a  three-dimensional  PML  in  spherical  coordinates  for  the  time-harmonic 
linear-acoustics  equations  according  to  Eq.  (2.26)  with  complex-coordinate  stretching  in  the  r  coor- 
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Figure  2:  Acoustic  scattering  of  a  plane  wave  on  a  unit  cylinder  with  a  cylindrical  PML.  Left:  Real 
component  of  the  numerical  solution.  Right:  Real  component  of  the  error  function.  Both  the  numerical 
solution  and  the  error  function  were  computed  on  the  hp  mesh  given  in  Fig.  1  (right). 


dinate  only.  To  assess  the  effectiveness  of  the  PML  combined  with  /ip-adaptivity,  we  compare  its 
performance  to  the  one  of  Infinite  Elements;  see  e.g.  Refs.  [4,  12]. 

Below,  we  present  numerical  results  for  the  scattering  of  a  plane  wave  on  a  sphere  with  dimen¬ 
sionless  radius  r  =  A  =  1,  where  A  denotes  the  dimensionless  wavelength.  We  set  the  dimensionless 
wavenumber  k  =  2tt  and  the  direction  of  the  incident  wave  v  =  (0,  0,-1).  The  computational  domain 
is  truncated  by  a  spherical  PML  for  2  <  r  <  3.  In  particular,  we  have  employed  coordinate  stretching 
according  to 


Jr  0  <  r  <  2 

I  r-iC{T-2f  r>2 


(2.32) 


in  order  to  have  z  G  C®  globally,  and  the  constant  C  =  —  log  e,  where  e  is  the  machine  precision. 
The  PML  is  truncated  by  a  homogeneous  Dirichlet  boundary  condition  p  =  0  at  r  =  3,  and  the  rigid 
scatterer  is  represented  by  the  Neumann  boundary  condition  dp/dn  =  —dp™’^/dn  on  F^v,  i-e.  at  r  =  1, 
where  is  the  incident  plane  wave. 

Figs.  5  (left)  and  (right)  show  cutaway  views  of  two  initial  meshes  used  for  comparison.  Mesh  (a) 
simply  captures  the  geometry  of  the  problem  with  the  PML  in  the  outer  layer  of  elements.  Mesh  (b) 
includes  one  /i-refinement  in  the  radial  direction.  In  Fig.  6  (right),  we  plot  the  percent  relative  error 
in  the  iL^-seminorm  (in  a  logarithmic  scale)  versus  the  total  number  of  degrees-of-freedom  N  (in  the 
algebraic  scale  for  uniform  p  =  2, . . .  ,9  and  both  meshes  (a)  and  (b).  The  error  is  evaluated 

only  over  the  domain  of  interest  1  <  r  <  2  with  the  exact  solution  (see  e.g.  [19,  sec.  2.1.2])  used  as  a 
reference.  We  see  that  the  radial  /i-refinement  is  essential,  and  that  the  curve  for  mesh  (b)  “looks  like” 
exponential  convergence  with  respect  to  p  in  this  range.  We  also  show  convergence  for  the  p-method 
when  the  PML  in  mesh  (a)  is  replaced  by  infinite  elements  on  the  sphere  r  =  2.  The  exact  solution 
is  then  analytic,  and  we  observe  exponential  convergence  in  p,  but  with  a  faster  rate  than  in  the  case 
with  the  PML.  Finally,  we  apply  our  algorithm  for  fully-automatic  /ip-adaptivity  (see  [21])  with  PML 
and  using  mesh  (a)  as  the  initial  coarse  grid.  We  plot  the  error  for  the  first  nine  coarse  grids  and 
corresponding  fine  grids  generated  by  the  algorithm.  The  final  coarse  grid,  shown  in  Fig.  6  (left), 
achieves  0.6%  relative  error  with  lOK  degrees-of-freedom,  and  the  corresponding  fine  grid  has  0.008% 
relative  error  with  137K  degrees-of-freedom.  We  note  that  the  algorithm  has  automatically  selected 
radial  /i-refinements  within  the  PML,  and  that  the  convergence  for  the  coarse  grid  appears  to  recover 
the  rate  delivered  by  infinite  elements.  In  Figs.  7  (left)  and  (right),  we  show  the  real  component  of 
the  solution  and  of  the  error  function,  respectively. 
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Figure  3:  Acoustic  scattering  of  a  plane  wave  on  a  unit  cylinder  with  a  Cartesian  PML.  Left:  hp- 
refined  mesh  corresponding  to  an  error  estimate  of  1%  and  consisting  of  4779  degrees-of-freedom. 
Right:  Real  component  of  the  numerical  solution  computed  on  the  hp  mesh  given  in  Fig.  3  (left). 

3.  Elasticity 

In  Section  3.1,  we  derive  the  PML  formulation  of  the  equations  governing  linear  elasticity.  In  Sec¬ 
tion  3.2,  we  present  numerical  results  that  demonstrate  the  validity  of  the  PML  formulation  for  a 
two-layered  elastic  medium  and  the  improvement  in  PML  performance  obtained  by  hp-adaptivity. 

3.1  PML  formulation  of  the  elasticity  equations 

In  this  section,  we  present  the  PML  formulation  for  linear  elasticity  and  discuss  its  implications.  To 
this  end,  let  us  recall  the  equations  of  linear  elasticity  in  K"  which,  in  time-harmonic  form,  are  given 

by 

{—puf^Ui  —  aijj  =  0  i  =  1, . . .  ,n 

1  ’  (3.33) 

^ij  C/ijkl  2  i.^k.1  L  j  1 ,  .  .  .  ,  71  , 

where  Ui  denotes  the  z-th  component  of  the  displacement  vector,  (Jij  is  the  stress  tensor,  p  is  the 
density  and  Eijki  is  the  elasticity  tensor  satisfying  the  usual  symmetry  properties 

Eijkl  —  Ejikh  E^jj^i  —  Eijlk-;  E^jkl  —  Ej^Hj  .  (3.34) 

The  second  symmetry  property  implies 

Eijki  ^{uk,i  +  upk)  =  EijkiUk,i  .  (3.35) 

For  an  isotropic  homogeneous  material,  the  elasticity  tensor  depends  on  two  constants  only,  viz. 

Eijki  —  p(.^ik^jl  4“  ^il^jk')  4“  ^^ij^kl  ;  (3.36) 

where  p  and  A  are  the  Lame  constants  and  5ij  is  the  Kronecker  delta.  In  Eqs.  (3.33)-(3.36)  and 
throughout  this  section,  we  make  use  of  the  Einstein  summation  convention. 

The  system  of  equations  (3.33)  can  be  complemented  by  various  boundary  conditions  of  which  we 
shall  restrict  ourselves  to  the  simplest  ones: 


Ui  =  0,  7  =  1,...,  77  onFu, 

^ij^j  —  ^  —  1,...,77  On  F  7\r  , 


(3.37a) 

(3.37b) 
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Figure  4:  Acoustic  scattering  of  a  plane  wave  on  a  unit  cylinder.  Left:  Real  part  of  the  error  estimate 
for  the  case  with  Cartesian  PML.  Right:  Percent  relative  error  versus  number  of  degrees-of- freedom 
for  the  cylindrical  PML  (exact  error  and  error  estimate)  and  for  the  Cartesian  PML  (error  estimate). 


where  rij  are  the  components  of  the  outward  normal  unit  vector,  and  gi  are  prescribed  tractions. 
Eq.  (3.37a)  prescribes  zero  displacements  and  corresponds  to  homogeneous  Dirichlet  boundary  con¬ 
ditions,  whereas  Eq.  (3.37b)  prescribes  given  tractions  and  corresponds  to  Neumann  boundary  condi¬ 
tions. 

To  derive  the  standard  variational  formulation  in  terms  of  the  displacement  vector,  we  multiply 
the  momentum  equations  (3.33) i  with  a  test  function  v  :=  Vi  G  V,  integrate  over  the  domain  LI  and, 
upon  integration- by-parts,  we  obtain 

/  dfl  —  cjijUjVi  dS  —  up'  /  pUiVi  dfl  =  0 ,  (3.38) 

where  we  have  tacitly  restricted  ourselves  to  test  functions  that  satisfy  homogeneous  boundary  con¬ 
ditions  on  the  Dirichlet  portion  of  the  boundary,  P Substitution  of  the  constitutive  law  (3.33)2  and 
Neumann  boundary  condition  (3.37b)  into  Eq.  (3.38)  then  yields  the  final  variational  formulation 

(  ugV 


'  /  EijkiUkjVij  dLt-  up  /  pUiVi  dLt=  giVi  dS,  (3.39) 

Vu  G  V , 

where  V  is  the  space  of  test  functions 

V  :=  {v  G  X  :  V  =  0  on  P^,}  ,  (3.40) 

which  constitutes  a  subspace  of  the  energy  space  X  := 

Complex-coordinate  stretching  of  the  linear-elasticity  equations  and  derivation  of  the  PML  for¬ 
mulation  is  straightforward,  and  we  shall  thus  discuss  only  the  Cartesian  form  here.  Under  analytic 
continuation  according  to  Eq.  (2.6)  with  xj  —>■  Zj{xj),  the  elasticity  equations  (3.33)  transform  into 


-PZ  LO  - =  0 


i  =  1, .  . .  ,n 
i,j  =  l,...,n. 


(3.41) 
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Figure  5:  Acoustic  scattering  of  a  plane  wave  on  a  sphere.  Two  different  initial  grids  for  assessing  the 
performance  of  a  spherical  PML  and  infinite  elements. 


where  Eq.  (3.41)i  has  been  multiplied  by  z’  :=  which  facilitates  the  subsequent  integration- 

by-parts,  since  the  prefactor  z' / z'^  is  independent  of  Xj. 

The  variational  formulation  of  (3.41)  can  then  be  derived  along  the  same  lines  as  for  the  standard 
case  considered  in  the  beginning  of  this  section,  which  yields 

'  MS  V, 

<  /  EijkiUk,iVij  dVl-  up  /  puiVi  dVt=l  giVi  dS,  (3.42) 

Vu  e  V , 

where  the  “modified”  elasticity  tensor  and  density  are  given  by 

Eijki  ■=  Eijki—j—^  (no  summation) 
p  :=  pz'  . 


(3.43a) 

(3.43b) 


In  (3.42),  the  test  space  V  is  a  subspace  of  the  energy  space  X  for  the  stretched  variational  form  (3.42) 
which  can  be  defined  as 

V  :=  {u  G  :  t;  =  0  on  F^,}  (3.44) 

with 

X=\v  :  /  \E^jki\vk.iVi,j  dfl  +  /  \p\v^v^dn  <  ooi  (3.45) 

I  JQ  JQ  J 

with  the  overbar  denoting  the  complex  conjugate.  Hence,  the  energy  space  (3.45)  corresponding  to 
the  variational  form  (3.42)  is  a  weighted  space  with  weights  that  incorporate  the  derivatives  of 
the  stretching  functions.  Note  that  under  complex-coordinate  stretching  the  minor  symmetries  of  the 
elasticity  tensor  (cf.  (3.34))  are  lost,  which  renders  the  elasticity  tensor  in  the  PML  non-physical. 
However,  since  the  major  symmetry  is  retained,  the  bilinear  form  (3.42)  is  complex-symmetric. 

Finally,  let  us  point  out  that  the  formulation  (3.42)  is  valid  also  for  layered  media  with  inter¬ 
faces  that  are  aligned  with  the  Cartesian  axes.  This  is  a  consequence  of  the  fact  that  the  solution 
remains  analytic  in  terms  of  the  coordinate  parallel  to  the  interface.  We  present  numerical  results 
demonstrating  the  validity  of  (3.42)  for  layered  media  in  Section  3.2. 
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Figure  6:  Acoustic  scattering  of  a  plane  wave  on  a  sphere.  Left:  hp-refmed  coarse  mesh  for  the 
spherical  PML  after  9  iterations  corresponding  to  0.6%  relative  error  and  consisting  of  lOK  degrees- 
of-freedom.  Right:  Comparison  of  /ip-coarse  and  fine  grids  with  the  p-method  for  spherical  PML  and 
with  the  p-method  for  infinite  elements. 


3.2  Numerical  results  for  an  elastic  layered  medium 

In  this  section,  we  demonstrate  the  validity  of  the  elasticity  PML  formulation  (3.42)  for  layered  media 
with  interfaces  that  are  aligned  with  the  Cartesian  axes.  To  this  end,  we  consider  a  vibrating  cylinder 
in  a  two-layered  medium.  Elastodynamic  wave  propagation  in  the  presence  of  material  layers  is 
especially  relevant  for  borehole  geophysics  applications;  see  e.g.  Ref.  [22]. 

We  solve  the  linear-elasticity  equations  on  a  rectangular  domain  (—4,4)^  with  a  circular  hole  of 
unit  radius  r  in  its  center.  The  domain  is  surrounded  by  a  Cartesian  PML  of  unit  thickness  in  which 
we  apply  complex-coordinate  stretching  according  to 

(  Xk  0  <  \xk\  <  4 

Note  that,  in  contrast  to  the  complex-coordinate  stretching  that  we  used  earlier  in  this  work  (see  for 
instance  Eq.  (2.31)),  the  present  form  of  stretching  accounts  also  for  the  term  aj{xj)  in  accordance 
with  Eq.  (2.6)  to  accelerate  the  decay  of  evanescent  waves  in  the  PML;  cf.  also  the  discussion  regarding 
the  complex  stretching  in  Section  2.2.  Such  evanescent  waves  are  generated  at  the  interface  between 
the  two  material  layers.  This  term  in  the  stretching  formula  warrants  that  evanescent  waves  are 
made  to  vanish  in  the  PML  before  reaching  the  domain  boundary.  We  remark  that  void  of  this  term, 
spurious  reflections  of  evanescent  waves  off  the  domain  boundary  would  compromise  convergence  and 
accuracy. 

Traction  boundary  conditions  corresponding  to  a  unit  pressure  are  applied  along  the  circumference 
of  the  hole,  and  homogeneous  Dirichlet  boundary  conditions  are  applied  on  the  outer  domain  boundary 
adjacent  to  the  PML.  The  material  data  are  given  in  terms  of  the  Poisson  ratio  v,  the  dimensionless 
density  p  and  the  dimensionless  Young’s  modulus  E 


V  =  0.3  ,  p  =  1 ,  E  = 


4  X2  >  0 

1  X2  <  0 


(3.47) 


Moreover,  we  set  w  =  tt. 

Starting  from  the  same  initial  mesh  as  displayed  in  Fig.  1  (left),  an  hp-refmed  mesh  correspond¬ 
ing  to  an  error  estimate  of  1%  in  the  iL^-seminorm  and  consisting  of  6508  degrees-of-freedom  is 
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Figure  7:  Acoustic  scattering  of  a  plane  wave  on  a  sphere.  Real  component  of  the  numerical  solution 
(left)  and  real  component  of  the  error  function  (right)  computed  on  the  hp  mesh  given  in  Fig.  6  (left) 
and  using  the  spherical  PML. 


obtained  after  16  iterations;  see  Fig.  8  (left).  Observe  the  refinements  that  are  selected  by  the  adap¬ 
tive  algorithm,  in  particular  in  the  PML  layer  to  resolve  the  enforced  rapid  decay  of  the  solution. 
In  Fig.  8  (right),  we  compare  the  convergence  of  /ip-adaptive  refinements  and  the  simpler  /i-adaptive 
refinements  with  approximation  degree  p  =  2  by  plotting  the  percent  relative  error  estimate  in  the  H^- 
seminorm  (in  a  logarithmic  scale)  versus  the  total  number  of  degrees-of-freedom  N  (in  the  algebraic 
scale  The  error  estimate  is  obtained  using  the  fine-grid  solution  as  a  reference  solution  and  is 

evaluated  over  the  entire  computational  domain  including  the  PML.  We  remark  that  the  /i-adaptive 
refinement  strategy  allows  for  anisotropic  refinements  that  are  crucial  for  resolving  the  PML-induced 
gradients  of  the  solution.  Fig.  8  (right)  indicates  that  both  hp-  and  /i-adaptive  refinements  deliver 
exponential  convergence.  Indeed,  in  the  preasymptotic  range,  /i-refinements  can  deliver  exponential 
convergence  before  algebraic  convergence  sets  in;  cf.  Ref.  [11,  ch.  14.3  and  14.4].  Note,  however, 
that  the  convergence  rates  and  error  levels  of  hp-  and  /i-refinements  are  significantly  different.  Hence, 
for  a  user-specified  error  tolerance,  an  /ip-adaptive  discretization  is  computationally  much  cheaper 
than  an  /i-adaptive  discretization.  Moreover,  it  is  obvious  that  for  the  resolution  of  PML-induced 
solution  gradients  an  h-uniform  discretization  would  be  prohibitive  in  terms  of  computational  cost. 
Our  numerical  results  thus  demonstrate  that  by  means  of  ft.p-adaptivity  reflections  from  the  PML  can 
be  minimized  to  an  arbitrary  tolerance  while  retaining  optimal  accuracy  for  the  least  computational 
expense.  In  Figs.  9  (left)  and  (right),  we  show  the  real  component  of  the  solution  of  the  horizontal 
and  vertical  displacement,  respectively,  obtained  on  the  /ip-mesh  given  in  Fig.  8  (left).  From  Fig.  9 
it  is  apparent  that,  in  comparison  with  the  greater  Young’s  modulus  for  X2  >  0,  the  smaller  Young’s 
modulus  for  X2  <  0  results  in  a  shorter  wavelength  of  the  solution  and,  thus,  necessitates  a  finer  hp 
discretization  (cf.  Fig.  8  (left)). 

4.  Maxwell’s  equations 

In  Section  4.1,  we  derive  the  PML  formulation  of  Maxwell’s  equations  in  two-dimensional  Carte¬ 
sian  coordinates.  In  Section  4.2,  we  present  numerical  results  that  demonstrate  the  improvement 
in  performance  of  the  PML  for  Maxwell’s  equations  that  can  be  obtained  by  using  an  /ip-adaptive 
discretization. 
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Figure  8:  Wave  propagation  in  an  elastic  two-layered  medium.  Left:  hp-mesh  corresponding  to  an 
error  estimate  of  1%  and  consisting  of  6508  degrees-of-freedom.  Right:  Percent  relative  error  estimate 
versus  number  of  degrees-of-freedom  for  /ip-adaptive  and  /i-adaptive  refinement. 


4-1  PML  formulation  of  Maxwell’s  equations  in  2D  Cartesian  eoordinates 

To  derive  the  PML  formulation  for  the  two-dimensional  Maxwell’s  equations,  let  us  recall  the  time- 
harmonic  form  of  Ampere’s  law  and  Faraday’s  law  in  three  dimensions 

{ioje  +  (j)E-V  xH  =  -  J™p  (4.48a) 

iujfiH  -I-  V  X  E  =  0  ,  (4.48b) 


respectively.  In  (4.48),  E  and  H  denote  the  electric  and  magnetic  field  vector,  respectively,  cu  is  the 
angular  frequency,  p,,  e  and  a  are  the  permeability,  permittivity  and  conductivity,  respectively, 
is  a  given  impressed  volume  current,  and  i  denotes  the  imaginary  unit. 

Eqs.  (4.48)  invoke  the  curl  operator  which,  in  Cartesian  coordinates,  is  expressed  as 


fdEs 

dE2 

dEi 

dEs 

dE2 

^E^\ 

\dx2 

dxs 

’  dX3 

dxi 

’  dxi 

dx2  ) 

(4.49) 


In  this  section,  we  shall  confine  ourselves  to  Maxwell’s  equations  in  two  dimensions.  However,  for 
future  reference,  we  include  the  PML  derivation  for  the  three-dimensional  Maxwell’s  equations  in 
Appendix  A.  For  a  vector- valued  two-dimensional  electric  field  E  =  {Ei,  E2,0)  with  components  that 
depend  only  on  xi,  X2,  only  the  third  component  of  the  three-dimensional  curl  operator  is  different 
from  zero 


V  X  E  ^ 


(4.50) 


Accordingly,  we  shall  define  a  two-dimensional  scalar-valued  curl  operator  as 


curl£l  := 


dE2 

dxi 


dEi 

dx2  ' 


(4.51) 


Similarly,  for  a  scalar-valued  one-dimensional  magnetic  field  H  =  (0,  0,  Hf)  the  curl  operator  reduces 
to 


V  X  H  = 


dx2  ’  dxi 


(4.52) 
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Figure  9:  Wave  propagation  in  an  elastic  two-layered  medium.  Real  component  of  the  horizontal 
displacement  (left)  and  of  the  vertical  displacement  (right). 


Invoking  Eqs.  (4.49)-(4.52),  the  two-dimensional  Maxwell’s  equations  in  time-harmonic  form  follow 
from  (4.48) 


itueEi  -I-  aEi  — 


iojeE2 


<jE2 

dE2 

dxi 


iufiHs  +  - - - —  =  0 


dHs 
dx2 
dHs 
dxi 
dEi 
dX2 


=  -  J^P 


=  -jr 


(4.53) 


If  the  solution  is  analytic  in  Xi  and  X2,  the  Cartesian  coordinates  Xi  and  X2  can  be  replaced  by 
the  complex  coordinates  zi  and  Z2^  respectively,  and  the  derivatives  can  be  interpreted  in  the  complex 


sense.  Under  complex-coordinate  stretching  according  to  (2.6)  with  Xj 
into 


Zj{xj),  Eqs.  (4.53)  transform 


(zwe  -I-  a)Ei  — 


iujg.H2, 


(iioe 
1 


a)E2 

5(4^2) 


dxi 


1  d{z[H3) 

z'yZ'2  dx2 

1  5(4^3) 

z(z2  dxi 
1  d{z[E^) 
z\  z(  dx2 


-J^p 


-J”P 


(4.54) 


=  0. 


where  we  have  made  use  of  the  fact  that  the  j-th  complex  coordinate  depends  only  on  the  j-th  real 

coordinate,  Zj  =  Zj{xj).  Upon  multiplication  of  Eq.  (4.54)i  by  iuJz[z2Fi  and  Eq.  (4.54)2  by 

with  F  :=  (El,  7^2)  a  test  function,  integration  over  domain  U  and  integration- by-parts,  we  obtain 


[  i-oj^e  +  ioja)  (^{z[Ei){z[Fi)  +  ^{z!2E2){z'2F2)]  dU  +  [  icoH^ 

Jn  V^i  ^2  /  Jn 


d{z[Fi)  d{z'2F2) 


dX2 


9a;  1 


dFt 


+  /  zwi73(^iEi(-n2)  +  4^’2ni)dE=-*oa  /  (z'j‘“P(z(Fi)  +  z'iJ™P(4E2))dU,  (4.55) 


Ida 


where  n  —  (711,712)  denotes  the  outward  normal  unit  vector.  Integrability  considerations  imply  that 
the  “generalized  tangential  component”  in  Eq.  (4.55), 


Ft  z'lFi{—n2)  +  ^2^2771 


(4.56) 
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must  be  continuous  across  inter-element  boundaries. 

Substituting  Eq.  (4.54)3  into  (4.55)  and  assuming  firstly,  a  nonhomogeneous  Dirichlet  boundary 
condition  on  the  physical  boundary 

nx  E  =  -nx  (4.57) 

with  £1“^^  a  prescribed  incident  wave,  secondly,  a  homogeneous  Dirichlet  boundary  condition  on  the 
PML  boundary 

z[Ei{-n2)  +  =  0  ,  (4.58) 

and  thirdly,  a  homogeneous  boundary  condition  for  the  test  function  on  the  entire  boundary 

Ft  =  0,  (4.59) 


we  obtain  the  final  variational  statement 


1  fdiz'^E^)  d{z[E,)\  (d{z'^F2)  d{z[F{) 


In 


dxi 


dxo 


dxi 


dx2 


(El 


-  [  {toh-iua)  (^iz[Ei){z[Fi)  +  ^{z'^E2)iz!,F2))  dn 

Jn  \^i  ^2  / 

=  —ito  f(  z' J™P(z(Fi)  -h  z'l  J™p(z'F2)  )  dfl  VF  .  (4.60) 


We  now  redefine  the  variables 


ZjEj,  j  =  l,2. 


(4.61a) 


z'F- 


J  =  l,2. 


(4.61b) 


Note  that  the  factor  z'  grows  algebraically  in  the  PML,  whereas  the  solution  Ej  decays  exponentially. 
Therefore,  the  redefined  variables  will  eventually  also  decay  exponentially  in  the  PML,  but  an  initial 
algebraic  growth  may  be  observed.  In  terms  of  the  redefined  variables,  Eq.  (4.60)  assumes  the  standard 
form 

E  G  E  +  E 


— \—r  cur  IF  curlF  —  (w^e  —  iua) 
fxz[z!2 


=  —lOJ 


(4.62) 


z^ +  z;4“PF2 j  dfl 
VEgE. 


Note  that  problem  (4.62)  is  complex  symmetric.  In  Eq.  (4.62),  E  denotes  a  lift  of  the  Dirichlet  data, 
and  the  test  space  F  is  a  subspace  of  the  “energy  space”  X  for  the  stretched  variational  form  (4.62) 
which  can  be  defined  as 


T  —  {E  gX  ■.  Ft  =  0  on  90} 


(4.63) 


with 


X 


1 


z  z: 


1^21 


curlF  G 


(4.64) 


Note  that  the  “energy  space”  AT  is  a  weighted  H (curl)  space  with  weights  associated  with  the  deriva¬ 
tives  of  the  stretching  functions. 

We  emphasize  that  the  redefinition  of  the  variables  according  to  (4.61)  is  not  merely  a  convenience. 
Integrability  assumptions  in  conformity  with  (4.64)  imply  continuity  of  the  tangential  component 
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Ft  of  the  redefined  variable  across  interelement  boundaries  and,  consequently,  the  redefinition  (4.61) 
enables  discretization  with  if  (curl)-conforming  elements.  Moreover,  let  us  point  out  that  the  redefined 
variables  satisfy  the  original  Maxwell’s  equations  and,  consequently,  the  variational  statement  (4.62) 
assumes  the  usual  form  associated  with  Maxwell’s  equations.  Therefore,  this  form  of  PML  is  typically 
also  referred  to  as  Maxwellian  PML;  see  e.g.  Ref.  [29]. 

4-2  Numerical  experiments 

In  this  section,  we  demonstrate  the  validity  of  the  PML  formulation  for  Maxwell’s  equations,  Eq.  (4.62), 
and  the  benefit  of  /ip-adaptivity  for  improving  the  PML  performance.  To  this  end,  we  consider  the 
scattering  of  a  plane  wave  on  a  cylinder. 

We  solve  Maxwell’s  equations  on  a  rectangular  domain  (—4, 4)^  with  a  circular  hole  of  unit  radius 
in  its  center.  The  domain  is  surrounded  by  a  PML  of  unit  thickness  in  which  complex-coordinate 
stretching  according  to 


Zk{Xk) 


Xk  0  <  \xk\  <  4 


(4.65) 


is  invoked.  Dirichlet  boundary  conditions  corresponding  to  an  incident  plane  wave 
are  applied  along  the  circumference  of  the  cylinder,  and  homogeneous  Dirichlet  boundary  conditions 
are  applied  on  the  outer  domain  boundary  adjacent  to  the  PML.  The  physical  parameters  are  given  in 
Table  1,  where  v  and  p  denote  the  direction  and  the  polarization  of  the  incident  wave,  respectively,  and 
r  is  the  cylinder  radius.  With  the  parameters  specified  in  Table  1,  the  speed  of  light  is  c  =  Ij ^/Jle  =  1 
and  the  wavenumber  is  fc  =  w/c  =  27r. 


Table  1:  Physical  parameters. 


CO  r  V  p  p,  e  a 

2tt  1  (-1,0)  (0,1)  1  1  0 


Starting  from  the  same  initial  mesh  as  displayed  in  Fig.  1  (left),  an  ft,p-refined  mesh  corresponding 
to  an  error  estimate  of  1%  in  the  ff(curl)-seminorm  and  consisting  of  14664  degrees-of-freedom  is 
obtained  after  16  iterations;  see  Fig.  10  (left).  Note  the  refinements  that  are  selected  by  the  adaptive 
algorithm,  in  particular  in  the  PML  layer.  In  Fig.  10  (right),  we  show  the  real  part  of  the  second 
component  of  the  solution  obtained  on  the  /ip-mesh  given  in  Fig.  10  (left).  The  corresponding  error 
estimate,  i.e.  the  difference  between  the  coarse  and  fine  grid  solution,  is  shown  in  Fig.  11  (left). 
We  remark  that,  in  contrast  to  acoustics  or  elasticity  problems,  the  solution  to  the  PML-modified 
Maxwell’s  equations  may  exhibit  an  initial  growth  in  the  PML  domain.  To  undermine  this  claim,  we 
display  in  Fig.  11  (right)  the  real  part  of  the  second  component  of  the  solution  for  a  wavenumber  k  =  1. 
The  initial  growth  in  the  PML  domain  is  caused  by  the  fact  that  we  have  chosen  to  discretize  the 
modified  variable  {z(Ei,  Z2E2)  rather  than  {Ei,E2)  directly.  In  the  PML,  the  stretching  coefficients 
Zi  initially  grow  algebraically,  whereas  the  solution  Ei  decays  exponentially.  Therefore,  the  modified 
variable  {z[Ei,  Z2E2)  also  grows  initially  before  it  exhibits  the  exponential  decay. 

5.  Conclusion 

In  this  work,  we  improved  the  performance  of  the  Perfectly  Matched  Layer  by  using  an  automatic  hp- 
adaptive  discretization.  The  PML  has  the  remarkable  property  of  having  a  zero  reflection  coefficient 
for  all  angles  of  incidence  and  all  frequencies  on  the  continuum  level.  However,  this  property  is 
compromised  under  discretization,  i.e.  at  a  discrete  PML  spurious  reflections  typically  do  occur.  Only 
for  vanishing  mesh  size,  i.e.  upon  convergence  to  the  continuum  solution,  the  property  of  having  a 


A.  PML  formulation  of  Maxwell’s  equations  in  3D  Cartesian  coordinates 


20 


p=l 


+1.34 

M 


H 

-1.32 


Figure  10:  Electromagnetic  scattering  of  a  plane  wave  on  a  cylinder.  Left:  Zip-refined  mesh  corre¬ 
sponding  to  an  error  estimate  of  1%  and  consisting  of  14664  degrees-of-freedom;  k  =  27r.  Right:  Real 
part  of  the  second  component  of  the  numerical  solution;  k  =  2tt. 


zero  reflection  coefficient  can  be  recovered.  By  means  of  an  automatic  Zip-adaptive  discretization,  we 
obtain  a  sequence  of  discrete  solutions  that  converges  exponentially  to  the  continuum  solution.  This 
allows  us  to  minimize  reflections  from  the  discrete  PML  to  an  arbitrary  level  of  accuracy.  The  strength 
of  our  Zip-adaptive  discretization  is  that  it  can  effectively  resolve  the  strong  PML-induced  gradients  of 
the  solution  by  adapting  the  element  size  Zi  and  polynomial  approximation  order  p  for  any  given  PML 
profile.  Since  the  Zip-adaptivity  is  automatic,  no  interaction  with  the  user  is  required.  This  renders 
tedious  parameter  tuning  of  the  PML  obsolete. 

We  demonstrated  the  improvement  of  the  PML  performance  by  Zip-adaptivity  and  the  versatility  of 
this  approach  by  numerical  results  for  acoustic,  elastodynamic  and  electromagnetic  wave-propagation 
problems  in  the  frequency  domain  and  in  different  systems  of  coordinates.  Our  results  show  that 
Zip-adaptivity  minimizes  reflections  from  the  PML  to  an  arbitrary  level  of  accuracy  while  retaining 
optimal  computational  efficiency.  Moreover,  we  showed  that  Zip-adaptivity  provides  the  freedom  to 
select  the  geometry  of  the  PML  independently  of  the  geometry  of  the  initial  mesh,  i.e.  the  respective 
geometries  need  not  coincide.  Our  results  suggest  that  Zip-adaptivity  is  ideally  suited  for  improving 
the  performance  of  the  PML  and,  conversely,  the  PML  is  the  ideal  absorbing  boundary  condition  for 
Zip-adaptive  finite-element  methods. 
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A.  PML  FORMULATION  OF  MAXWELL’S  EQUATIONS  IN  3D  CARTESIAN  COORDINATES 
As  an  extension  of  Section  4.1,  we  derive  in  the  sequel  the  PML  formulation  for  Maxwell’s  equations 
in  three  dimensions  and  Cartesian  coordinates.  To  this  end,  it  is  convenient  to  rewrite  the  time- 
harmonic  form  of  Ampere’s  and  Faraday’s  law,  Eqs.  (4.48a)  and  (4.48b),  respectively,  expressing  the 
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Figure  11:  Electromagnetic  scattering  of  a  plane  wave  on  a  cylinder.  Left:  Real  part  of  the  second 
component  of  the  error  estimate;  k  =  2tt.  Right:  Real  part  of  the  second  component  of  the  numerical 
solution;  k  =  1. 


curl  operator  by  means  of  the  permutation  symbol  Sijk  that  is  defined  as 

{0  ,  if  any  two  indices  are  repeated 

1 ,  if  ijk  is  an  even  permutation  of  123  , 

—1 ,  if  ijk  is  an  odd  permutation  of  123 

which  then  yields 

{iuJe  +  a)E,-J2s^JkHk,J  =  -  4“^  f  =  l,2,3 
^  ^  ^kmn^n,m  =0  =  1,  2,  3  , 

m,n 


(A.66) 


(A.67) 


respectively.  Mind  that  in  this  section,  we  write  sums  explicitly  using  the  summation  symbol,  rather 
than  making  use  of  the  Einstein  summation  convention.  Under  complex-coordinate  stretching  accord¬ 
ing  to  Eq.  (2.6)  with  Xi  Zi{xi),  system  (A.67)  transforms  into 

(iuje  +  a)Ei  -  ^ (A. 68a) 
iUJfJjH]^  +  ^  ^  E/ri,m  —  0  ■  (A. 68b) 

m,n 

Upon  multiplication  of  Eq.  (A. 68a)  by  z'  :=  we  obtain 


{iuje  +  a)^{ZiEi)  =  -z'Jf^'^ 

*  3,k 


(A.69) 
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We  then  multiply  Eq.  (A. 69)  by  a  test  function  Fi,  integrate  over  a  domain  Cl  and  integrate  by  parts, 
which  yields 


f  i{iuje  +  a)-^{ZiE,){z'iF,)+'^e^jk^Hk{zlF,)A  dn 
Jn  y  ^  J 

=  -  [  [  y^£ijk^Hk{zlF,)njdS.  (A.70) 

Jn  Jan  -  j. 

Next,  we  multiply  Eq.  (A.70)  by  iw  and  eliminate  Flk  from  the  volume  integral  by  invoking  (A. 68b) 

r  j  ^  ^  \ 

/  I  e  +  iwcr)  ^  {z^Ei){Zj^Fi)  ,  ,  (ZndJn)  ,m{z^Ei)  J  j  dfl 

Jn  y  y^i)  d-  -  ^  J 

=  -luj  [  ^ji^^iz'Ei)  dn  +  luj  f  V e.jk^jHkiz'F^hj  dS  .  (A.71) 
Jn  JdQ  z^z^ 

Summation  over  index  i  and  some  straightforward  algebraic  manipulations  then  yield 
/  I  (— +  iwcr)  ^  7^(2:jEj)(2:jEi)  —  —  ^  ekmn  £kij 

^  \  i  ^  ^  k  \m,n  n  m  /  \  ^,3  ^  ^  / 

=  -i^  !  y\  ^Jr"^i<F,)  dn  +  iiv  f  V  £^jk^Hk{z[F,)nj  dS  .  (A.72) 

Jn  Jan ■  j,  z^z. 


dO. 


Upon  redefining  the  variables 


Ei  <- 

-  z[E„ 

f  =  1,2,3, 

(A.73a) 

77,4- 

-  z'H,, 

f  =  1,2,3, 

(A.73b) 

-  z-F„ 

f  =  1,2,3, 

(A.73c) 

and  incorporating  into  Eq.  (A.72)  a  Neumann  boundary  condition  with  an  impressed  surface 
current,  we  obtain  the  final  variational  formulation 

f  E  eE  +  T 


a;^e  +  za;a)^  X]  (  X]  y  ^n,m  )  ( 


=  —lU) 


Z‘ 

z[ 


Eur-c" 


■f  iui 


{  yFeE, 

(A.74) 

where  E  denotes  a  lift  of  the  Dirichlet  data,  and  the  test  space  :F  is  a  subspace  of  the  “energy  space”  X 
for  the  stretched  variational  form  (A.74)  which  can  be  defined  as 


T  ■.=  {F  ^  X  :  rixF^OonEi)} 


(A.75) 
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Note  that  the  “energy  space”  X  is  a  weighted  H (curl)  space  with  weights  associated  with  the  deriva¬ 
tives  of  the  stretching  functions. 

Since  the  electric  and  magnetic  field  redefined  according  to  (A. 73)  satisfy  the  original  Maxwell’s 
equations  in  a  medium  with  complex  permittivity  and  permeability  constitutive  tensors,  this  form  of 
PML  is  typically  also  referred  to  as  Maxwellian  PML;  see  e.g.  Ref.  [29].  With  complex  permittivity 
and  permeability  constitutive  tensors  redefined  according  to  e  =  eA  and  jl  =  fxA  with 


A 


eiOi 


+  6262 


+  6363 


(A.77) 


such  PML  formulation  yields  a  physically  meaningful  material.  In  contrast,  the  PML  formulation  of 
the  linear-elasticity  equations  yields  a  modified  elasticity  tensor  that  is  non-physical;  cf.  Section  3.1. 
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